Heat conduction in a one-dimensional gas of elastically colliding particles of unequal 

masses 
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We study the nonequlibrium state of heat conduction in a one-dimensional system of hard point particles of unequal masses 
interacting through elastic collisions. A BBGKY-type formulation is presented and some exact results are obtained from it. 
Extensive numerical simulations for the two-mass problem indicate that even for arbitrarily small mass differences, a nontrivial 
steady state is obtained. This state exhibits local thermal equilibrium and has a temperature profile in accordance with the 
predictions of kinetic theory. The temperature jumps typically seen in such studies are shown to be finite-size effects. The 
thermal conductivity appears to have a very slow divergence with system size, different from that seen in most other systems. 

PACS numbers: 44.10.+i, 05.45.-a, 05.60.-k, 05.70.Ln 



Introduction: The problem of finding a one- 
dimensional system of interacting particles, evolving 
through Newtonian dynamics, in which Fourier's law of 
heat conduction holds, has attracted considerable inter- 
est in recent years Q . One of the first studies in the field 
was the work of Rieder et al || who obtained the exact 
nonequilibrium steady state (NESS) of a chain of cou- 
pled harmonic oscillators. They found a constant tem- 
perature profile in the system and a heat current that 
was independent of system size. This result was not 
too surprising, considering that the harmonic chain is 
integrable and there is no scattering of phonons. Since 
then, a large number of studies have looked at the ef- 
fects of introducing impurities, nonlinearity and external 
potentials [3 — 8]. The specific question asked in most 
of these studies is the following. Defining the thermal 
conductivity K of a system through the linear response 
formula j = —K(T)WT, where j is the the current and T 
the local temperature, is there a one-dimensional model 
which gives a finite K ? Contrary to initial expectations, 
it was found that adding impurities and/or nonlinearity 
into a system did not result in finite K. The first model 
in which a finite K was found was the so-called ding-a- 
ling model Q in which alternate particles on the line are 
bound to harmonic springs and the other particles move 
freely between pairs of the bound particles. Numerical 
studies of several other models have also given finite con- 
ductivity ||. A common feature in these models is the 
presence of external potentials which lead to momentum 
nonconservation. Recently it has been proved rigorously 
|Tof that the conductivity as given by the Green-Kubo 
formula always diverges in one-dimensional momentum 
conserving systems with finite pressure. Two recent stud- 
ies [ pd| have reported finite conductivity in momentum- 
conserving systems but both of these have vanishing pres- 
sure and so there is no contradiction. 

While the proof of anamalous conductivity in momen- 
tum conserving systems has been a significant progress, 



there remain many issues that still need to be addressed. 
First of all, the proof uses the Kubo formula and this 
is not fully satisfactory since the validity of the Kubo 
formula, even in the limit of arbitrarily small tempera- 
ture gradients, has never been established rigorously. It 
may be noted that any derivation of the Kubo formula 
for thermal conductivity essentially makes the following 
assumptions: (i) the NESS is in local thermal equilib- 
rium (LTE) (ii) that Fourier's law is valid and (iii) re- 
gression of equilibrium fluctuations occur in the same 
way as nonequilibrium processes fl2|| . These assump- 
tions are physically motivated but may not be true in 
all cases. Secondly we note that while the focus of most 
of the work on heat conduction has been addressed to 
obtaining Fourier's law, the more general problem is one 
of understanding the nonequilibrium energy current car- 
rying state of a many-body system. For example one 
important question is the existence or otherwise of LTE 
in the steady state. Thus one would like to know if it 
is possible to define a slowly varying temperature field 
which determines all other local properties in the sys- 
tem. This point has not attracted much attention, even 
though it is quite crucial even for stating Fourier's law 
and using results of linear response theory. Naively one 
might expect that if thermal currents in the system van- 
ish in the thermodynamic limit, then LTE should hold 
but this has been shown not to be true always [Q. An- 
other interesting question is the determination of the 
temperature profile itself. These questions are clearly 
of interest and can be asked independently of whether or 
not Fourier's law holds. Finally we note that numeri- 
cal studies of heat conduction are problematic for several 
reasons. One needs long time numerical solution of non- 
linear differential equations which is time-consuming and 
not very accurate. This, in addition to long equilibration 
times typically occuring in such systems, restricts one to 
small system sizes. Also the treatment of boundary re- 
lated problems, such as that of temperature jumps, is 
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not straightforward. Hence it has often been difficult to 
arrive at correct conclusions from numerical studies and 
it is desirable to have more accurate studies. 

In this paper we study heat conduction in a system 
of hard elastic particles of unequal masses moving on a 
one-dimensional line. The only interaction between the 
particles is through elastic collisions. This model was 
first considered by Casati |[4| as a possible candidate for 
obtaining Fourier's law but the numerical results were 
insufficient to draw any definite conclusions. More re- 
cently this model has been studied by Hatano H who 
obtained a diverging conductivity. The equal-mass case 
with dissipation was studied by Du et al |^| who obtained 
a rather surprising NESS which implied a breakdown of 
usual hydrodynamics. 

Here we present extensive and accurate numerical work 
on this model and also an analytic formulation of the 
BBGKY-type. Our aim has been to give a more de- 
tailed characterization of the NESS than has been pre- 
viosly done. The present model is particularly suitable 
for this purpose for two reasons: (i) simulations of this 
model do not require a numerical solution of nonlinear 
differential equations and it is possible to obtain very 
accurate results, (ii) analytically the BBGKY hierarchy 
has a comparatively simpler structure and some exact 
statements can be made. The most interesting result 
obtained is that the steady states for the case of equal 
masses and the case with arbitrarily small mass differ- 
ences are completely different. In the thermodynamic 
limit the latter case exhibits LTE and the temperature 
profile approaches a form predicted by kinetic theory for 
a one-dimensional gas. This is surprising since kinetic the- 
ory predicts a finite conductivity with a T 1 / 2 dependence 
on temperature. On the other hand our model is momen- 
tum conserving and the proof for diverging Kubo con- 
ductivity holds. In our finite size studies we find a slow 
divergence of the conductivity (~ L a with a < 0.2). Our 
work also clarifies some of the problems related to bound- 
ary effects. The jumps in temperature at the boundaries 
are shown to be finite-size effects which are studied sys- 
tematically. 

Definition of model: We consider N point particles 
numbered i = 1,...N moving in a one-dimensional box 
extending from to L. The mass, position and velocity 
of the ith particle are denoted by m^, Xi and Ui. The 
only interaction between particles is through elastic col- 
lisions. After a collision between particle i and (i + 1), the 
new velocities are obtained from momentum and energy 
conservation and given by the linear equations: 



m l - m t+ i 
mi + m l+ i 
2m; 
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Ui+l 
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-u l+1 . (1) 

rn i+ i m l + m l+1 

Between collisions the particles travel with constant ve- 
locity. The coupling to heat baths is implemented by 



using Maxwell boundary conditions. Thus whenever a 
particle of mass m collides with a wall at temperature 
T, it is reflected back with a velocity chosen from the 
distribution P(u) = {m\u \ /T)exp(-mu 2 / (2T)). The 
temperatures at the two ends are taken to be Ti and T 2 . 

The case when all particles have equal masses can be 
solved easily and behaves similarly to the ordered har- 
monic springs cas e |^]. The "temperature" profile is fiat 
(with the value yTi^j, current is independent of sys- 
tem size and there is no LTE. The equal-mass problem 
is integrable and essentially reduces to a single-particle 
problem and so the results are not surprising. As soon 
as the masses are made different, the system becomes 
nonintegrable and is expected to have good ergodicity 
properties Jll| , and correspondingly the NESS should be 
very different. 

BBGKY-type equations: In principle, a complete so- 
lution of the heat conduction problem could be ob- 
tained from the steady-state solution of the master equa- 
tion for evolution of the ./V-point distribution equation 
p({xi,ui},t). In practice this is difficult and a sim- 
pler approach is to work with the so-called BBGKY 
hierarchy which deals with reduced distribution func- 
tions [ p^| . Let us first make the following definitions: 
pi(x,u,t) = (5(x - xi)d(u - ui)), pi ! i + i(x 1 ,u 1 ;x 2 ,u 2 ) = 
(6(x! - xi)5(ui - ui)8(x 2 - xi +1 )8{u 2 - u l+1 )), 
where (A) = J p({xi,ui},t)A({xi,ui})Y[idxidui. Fur- 
ther we define pi t i+\ (x, u\,u 2 ) as the number of collisions 
per unit time occuring at x between the Ith and (I + l)th 
particles with respective velocities U\ and u 2 . Clearly 
Pi,i+i{x,Ui,u 2 ) = pi ; i + i{x,u x ;x,u 2 ){ux - u 2 )6(u x - u 2 ). 
In terms of these the BBGKY-type equations relating 
one-point functions to two-point functions are the fol- 
lowing: 

dpi (x, u,t)/dt + udpi (a;, u, t)/dx = 

ui,u 2 )8(u' 2 - u)duidu 2 - / pi-ij(x,ui,u)dui 



Pl,l+\[x,Ui,U2]5{u' 1 - u)duidu 2 - p u+ i(x,u,u 2 )du 2 . (2) 



These equations hold for x ^ 0, L. At the boundaries, 
the distribution functions satisfy appropriate boundary 
conditions. 

The physical observables that we will be interested in 
are the particle density n{x,t), the energy density e{x,t) 
and the energy current density j{x,t). These can be 
expressed in terms of the one-point functions pi(x,u,t). 
Thus we have: 

n(x, t) — (S~2 8(x — xi)) = / pi(x 7 u,t)du 

. . ,\r-^ miuf yr-^ m l f 2 / \ 

I I J 

j( x ,t)={J2^j L S(x-x l )) = J2 1 Y J u 3 pi(x, Ul t)du (3) 
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The temperature field T(x) is defined as T(x, t) = 
2e(x, t) /n{x, t). Our first result is the following exact cur- 
rent conservation equation: de(x,t)/dt+dj(x,t)/dx = 0. 
This result is easily obtained by multiplying Eq. (||) 
by m;u 2 /2, integrating over u, and summing over all 
I. There is a pairwise cancellation of all terms on the 
right-hand side. Similarly by multiplying Eq. (|^) by 
to/w/2, integrating and summing gives, in the steady 
state: de(x)/dx — 0. Thus for any choice of masses 
{mi}, the energy density in the steady state is constant in 
space. Physically the constancy of energy density follows 
from the constancy of pressure and the linear relation 
between the two quantites in an ideal gas. However this 
does not imply a constant temperature profile since the 
temperature also depends on the number density n(x) 
which is not constant. 

From the fact that the dynamics [Eq. ([!])] is invariant 
under a constant scaling of the masses, it follows that the 
temperature profile does not change under mi — > vrrii. 
Also from the boundary conditions it is easily shown that 
T(yTi,vT2,x) — vT(Ti,T2,x). From now on, we shall 
consider the dimer case where we consider only particles 
of two different masses mi and mi placed alternately on 
the line. Because of the above scaling properties the only 
independent variables are the ratios, milm\ and T2/T1, 
and N. We will henceforth consider the case mi = 1, 
777,2 = to, T\ = 2 and T2 = 8. 

Numerical results: In our numerical simulations we let 
the system evolve with the appropriate boundary condi- 
tions and compute time-averages of various quantities in 
the steady state. The time evolution does not require nu- 
merical solution of differential equations since the exact 
solution is essentially known. The system is evolved by 
computing successive collision times and updating veloc- 
ities using the collision rules Eq. (^) . The only errors are 
those due to round-off. We have verified that the simula- 
tions reproduce all the exact known results both for the 
equal and the unequal mass cases. 

In our simulations we vary 5 = to — 1 and N. The 
mass values 8 = 0.078,0.11,0.22,0.44 were studied for 
particle numbers N = 41, 81, 161, 321, 641 and 1281. The 
size of the boxes were changed so that in all cases the 
average density of particles was fixed at 2. The number 
of particles is chosen to be odd so that at both ends the 
particles in contact with the bath have the same mass. 
The number of collisions over which the averaging is done 
was between 10 9 — 10 10 . In all cases we checked that 
increasing the time of averaging by a factor of 10 did not 
significantly change the data. 

In Fig. ([!]), we plot the steady state temperature pro- 
files at different system sizes for 8 = 0.22. The tem- 
perature has a smooth and continuously varying profile. 
There is a jump at the boundaries which decreases as the 
system size is increased and is expected to vanish in the 
thermodynamic limit when the current also becomes van- 
ishingly small. We find that this is true for any non-zero 



8. For smaller 8 one needs to go to larger system sizes 
to get the same temperature profile. Infact for small <5 
and large N the temperature-profile depends on 8 and N 
only through the scaling combination 8 2 N. This is illus- 
trated in Fig. (H) where we plot data corresponding to 
five different values of 8, each with a different N, chosen 
such that 8 2 N is the same. Note that for 5 = 0, the tem- 
perature profile (which is flat and given by T = \JT\T2) 
and the energy current are both independent of system 
size; thus 8 = is a singular point, while the steady state 
in the limit 8 — > 0, N — > 00 with 8 2 N constant is quite 
different and nontrivial. 

For fixed value of 8, as we increase N the tempera- 
ture profile approaches a limiting form. Quite amazingly 
this limiting form seems to be exactly one that would be 
predicted by kinetic theory. We recall that kinetic the- 
ory for a one-dimensional gas predicts Fourier behaviour 
with K ~ T 1 / 2 and hence a nonlinear temperature pro- 
file T k (x) = [T^ /2 {l^x/L)+T^ /2 x/L} 2 / 3 . This has been 
plotted in Fig. (|l|). We find that the following scaling 
form gives a reasonable collapse of our data: 

T(x,N,6)=T k (x) + -^pg(x). (4) 

The inset in Fig. ([!]) shows the collapse of data for 8 = 
0.22 obtained using the above scaling form with 7 = 0.67. 

We now look at the dependence of K on system size. 
In Fig. (|J) we plot j versus L for 8 = 0.22. It is clear from 
the data that the conductivity which is proportional to 
Lj has a slow divergence given by K ~ L a with a < 0.2. 
We note that this is significantly different from the sys- 
tem size dependence of K ~ L aA found by Hatano in the 
same model and also in other momentum-conserving 
systems like FPU and the diatomic Toda I^Q. 

Finally we have checked for LTE: to do so we compute 
the steady state expectation u^ 4 \x) = (J2i m i u t^( x ~ 
x{)). If there was LTE, this quantity would be deter- 
mined by the local temperature T(x). In Fig. (^) we plot 
w( 4 )(x), as determined directly by taking time averages 
and also the value predicted from the local temperature 
T(x). We see that at large system sizes LTE is indeed 
achieved. We also find that for smaller values of 8, one 
needs to go to larger system sizes to get LTE. 

In summary we have studied heat conduction in the 
unequal mass problem which appears to be the simplest 
nontrivial deterministic system, in one dimension, for 
which a very detailed investigation of the NESS can be 
made. Our study shows that a meaningful hydrodynamic 
description of the steady state is possible even in a situa- 
tion where (presumably) K — ► 00 in the thermodynamic 
limit. It is clear that further studies of this model can 
throw much light on the difficult problem of transition 
from the microscopic to macroscopic description in the 
context of nonequilibrium phenomena. 
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FIG. 1. Temperature profiles for 5 — 0.22 are plotted for 
system sizes N = 21(indicated by o), 41, 81, 161, 321, 641 and 
1281 (+). The solid line is the prediction of kinetic theory. 
In the inset we have plotted T sc = N 0S7 (T(x) - T k (x) [see 
Eq. (Eh] with the data for TV = 161, 321, 641 and 1281. 
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FIG. 2. Temperature profiles obtained with five different 
sets of values for 5 and TV with 5 2 N constant. 
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FIG. 3. Plot of j versus L for S = 0.22 The straight line 
shown corresponds to the decay j ~ 1/L 83 . 
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FIG. 4. Plot of u' ' as determined from a direct time 
averaging (NeqAv) and from the local temperature (EqAv) 
for two system sizes. For the bigger system, the curves almost 
coincide. 
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